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ABSTRACT 

The  returns  from  shallowly  buried  targets  measured  using  Ground  Pen¬ 
etrating  Radar  (GPR)  are  typically  obscured  by  a  strong  background  signal 
comprised  of  the  reflections  from  the  air-soil  interface.  A  Kalman  filter-based 
approach  is  proposed  to  estimate  this  background  signal  and  to  separate  it 
from  the  target  return.  In  the  absence  of  the  target  the  filter  operates  using 
a  “quiescent  state  model”  in  which  it  computes  the  background  estimate.  A 
statistic  based  on  measurement  innovation  is  applied  to  detect  the  target  posi¬ 
tion.  Upon  detection  the  state  is  augmented  by  a  new  component  which  allows 
for  the  change  of  the  signal  corresponding  to  the  presence  of  the  target  return. 
The  augmented  state  model  is  used  until  it  is  reverted  to  the  quiescent  model 
by  another  decision. 
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A  Kalman  Filter-Based  Approach  to  Target  Detection  and 
Target-Background  Separation  in  Ground  Penetrating 

Radar  Data 


EXECUTIVE  SUMMARY 

Ground  Penetrating  Radar  (GPR)  has  shown  promise  for  detecting  landmines  with 
minimal  or  no  metal  content.  A  problem  related  to  this  technology  is  that  the  returns 
from  shallow  buried  objects  can  be  obscured  by  the  ground  return.  In  the  scope  of  its 
countermine  project  DSTO  is  investigating  signal  processing  techniques  for  improving 
detection  of  shallow-buried  targets  within  the  ground  return. 

This  report  describes  a  Kalman  filter-based  approach  to  target  detection  and  clutter 
suppression  in  near  surface  GPR  data.  The  detection  problem  is  posed  as  the  one  of 
detecting  local  anomalies  in  the  soil  dielectric  half-space,  where  the  soil  properties  are 
assumed  to  slowly  change  across  the  half-space  as  a  function  of  distance.  Background 
estimation,  target  detection  and  target-background  separation  are  treated  as  mutually 
related  processes  and  are  performed  within  an  integral  Kalman  filter-based  computational 
procedure.  In  the  absence  of  the  target  the  filter  operates  using  a  quiescent  signal  model 
in  which  it  computes  the  background  estimate.  A  statistic  based  on  the  measurement 
innovation  obtained  from  the  Kalman  filter  is  applied  to  detect  the  target  position.  Upon 
the  detection  the  state  is  augmented  by  new  components  which  allow  for  the  change  in 
the  signal  corresponding  to  the  target  return. 

This  technique  was  tested  using  a  number  of  targets  measured  in  several  soil  environ¬ 
ments  and  proved  to  be  superior  to  the  standard  background  subtraction  method.  Further 
tests  are  required  in  order  to  determine  the  effectiveness  of  this  technique  in  the  environ¬ 
ments  that  contain  a  large  number  of  false  targets  (clutter).  Several  extensions  to  the 
current  method  which  are  expected  to  provide  an  increase  in  the  detection  capability  have 
been  suggested. 
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1  Introduction 


Impulse  time-domain  Ground  Penetrating  Radar  (GPR)  is  capable  of  detecting  buried 
targets  with  little  or  no  metal  content  and,  as  such,  has  emerged  as  a  promising  tool  for 
landmine  detection.  GPR  radiates  short-duration  pulses  of  electromagnetic  energy  into  the 
ground  and  records  backscattered  signatures  composed  of  the  reflections  from  the  target 
dielectric  surfaces  ( e.g .,  top  and  bottom  of  a  mine  plastic  casing)  and  internal  metallic 
components  (e.g.,  firing  pin).  However,  GPR  signatures  of  shallowly  buried  targets  such 
as  landmines  are  normally  obscured  by  a  background  signal  (clutter)  comprised  of  the 
reflections  from  the  ground  surface  and  the  antenna  crosstalk,  which,  is  some  cases,  may 
impede  detection. 

As  the  target  and  the  ground  returns  have  similar  spectral  characteristics  standard 
radar  clutter  suppression  techniques  can  not  be  used  for  their  separation  [5].  A  simple 
alternative  is  to  apply  background  subtraction  whereby  the  background  signal  is  estimated 
as  the  mean  of  the  unprocessed  ensemble  of  GPR  signals  using  either  all  signals  in  the 
ensemble  or  a  spatial  moving  average  filter  to  obtain  a  locally  adaptive  background  esti¬ 
mate  [4],  [5].  This  estimate  is  next  subtracted  from  the  unprocessed  GPR  signals  in  the 
ensemble.  The  resulting  data  is  usually  utilised  to  locate  the  position  of  a  possible  target 
which  may  be  followed  by  target  recognition  algorithms  [2].  We  note  that  the  techniques 
described  in  [4]  and  [5]  do  not  take  into  account  the  target  position  and  can  use  target 
return  signals  to  compute  background  estimate,  which,  in  some  cases,  may  modify  the 
distribution  of  the  target  signatures  after  the  subtraction.  By  contrast,  in  this  report 
background  estimation,  target  detection  and  target-background  separation  are  treated  as 
mutually  interrelated  processes  and  are  performed  withiri  an  integral  Kalman  filter-based 
computational  procedure.  A  Kalman  filter  approach,  motivated  by  the  variable  state  di¬ 
mension  (VSD)  method  described  in  [1],  is  utilized  to  estimate  the  background  signal  and 
to  separate  it  from  the  target  signal.  In  the  absence  of  the  target  the  filter  operates  using  a 
quiescent  state  model  in  which  it  computes  the  background  estimate.  The  target  is  defined 
as  a  local  anomaly  in  the  soil  dielectric  and  is  detected  as  an  abrupt  departure  from  the 
estimated  background  signal.  A  statistic  based  on  the  measurement  innovation  obtained 
by  the  Kalman  filter  is  applied  to  detect  the  target  position.  Upon  detection  the  state  is 
augmented  by  new  components  which  allow  for  the  change  of  the  signal  corresponding  to 
the  target  return.  The  augmented  state  model  is  used  until  it  is  reverted  to  the  quiescent 
model  by  another  decision. 

It  is  assumed  that  the  background  signal  is  a  slowly  changing  function  of  the  spatial 
position.  The  magnitude  of  this  signal  at  one  particular  depth,  and  for  incrementally 
updated  spatial  positions,  is  modelled  using  a  random  walk  model.  Based  on  this  model, 
the  noisy  measurements  of  the  background  signal  are  represented  using  a  state  space  form 
and  Kalman  filter  is  applied  to  obtain  the  recursive  state  estimates.  In  particular,  a  set 
of  Kalman  filters  is  run  on  horizontal  strips  of  GPR  data  which  are  adjacent  in  depth.  At 
the  spatial  positions  where  the  target  is  detected  the  state  is  augmented  by  the  new  state 
component  modelled  as  stochastic  bias.  This  component  is  used  to  compute  the  estimate 
of  the  unknown  target  signal  which  is  assumed  to  be  superimposed  onto  the  background 
signal  at  the  position  of  the  target. 

The  report  is  organized  as  follows.  Section  2  describes  the  Kalman  filter-based  ap- 
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proach  to  background  estimation  and  target-background  separation.  Section  3  then  presents 
our  approach  to  target  detection  using  measurement  innovation  and  Section  4  explains  the 
technique  that  enables  adaptation  of  the  algorithm  to  slow  changes  of  the  environment. 
The  experimental  results  are  presented  in  Section  5. 


2  Background  Estimation  and 
Target-Background  Separation  Using  Kalman 

Filter 

In  the  process  of  acquiring  GPR  signatures  the  transmit  and  receive  antennas  are 
moved  over  the  ground  surface  at  approximately  constant  velocity  v  and  height  h  so  as  to 
cross  the  buried  target.  The  received  backscattered  signals  (or  traces)  u(n,  k )  are  collected 
in  fixed  time  intervals  A tc  where  n  —  0, . . . ,  N  —  1  denotes  the  signal  time  samples  and 
k  =  0, 1, . . . ,  corresponds  to  the  position  of  the  receive  antenna,  vkAtc.  Using  the  additive 
signal  model  we  can  write: 

u(n,  k)  =  s4(n,  k)  +  sb(n,  k)  +  w(n,  k),  n  =  0,  —  1  (1) 

where  ^(n,  k)  is  the  target  signal,  sb(n,k )  is  the  background  signal  (superscript  t  stands 
for  “target”  and  b  denotes  “background”),  and  w(n,k)  is  additive  noise.  If  the  signal 
u(n,  k)  is  measured  over  the  area  with  no  target  present  ^(n,  A:)  =  0,  n  =  0, . . . ,  N  —  1. 

Two-dimensional  GPR  plots  are  obtained  by  stacking  the  traces  u(n,k )  measured  at 
consecutive  spatial  positions.  The  horizontal  axes  of  a  plot  corresponds  to  the  antenna 
position  and  the  vertical  axis  denotes  time  duration  of  the  backscattered  signal  (equiv¬ 
alently,  proportional  to  the  depth).  The  plots  are  composed  of  sections  of  traces  that 
contain  target  signal  and  those  measured  over  the  areas  with  no  target  present,  which 
constitute  the  background  signals  and  the  goal  is  to  detect  target  and  segment  the  data 
into  the  target  and  the  no-target  areas.  To  enable  efficient  processing,  GPR  plots  are  di¬ 
vided  into  non-overlapping  horizontal  strips  which  correspond  to  the  layers  of  the  ground 
that  are  approximately  constant  in  depth.  Each  strip  is  processed  separately  as  explained 
bellow. 

We  define  the  measurement  vector  as  u p(k)  =  [u(pm,  k),u(pm  + 1,  k), . . . ,  u((p+  l)m  — 
1,  k)]T,  where  p  =  0, 1, . . . ,  P  -  1,  and  P  =  [N/m]  {{z\  represents  the  largest  integer  <  z). 
That  is,  up(k )  equals  the  section  of  length  m  of  the  received  GPR  signal  u(n,  k )  which  starts 
from  the  pm-th  signal  sample  and  up(k),  k  =  0, 1, . . . ,  is  a  horizontal  strip  of  GPR  data 
of  width  m  positioned  at  pm.  Following  the  random  walk  model  for  the  background-only 
signals,  we  define  the  quiescent  state  using  the  following  linear  difference  state  equation: 

s bp(k)  =  sbp(k  ~  l)  +  (2) 

for  p  =  0, 1, . . . ,  P  —  1,  where  s b(k)  =  [ sb(pm ,  k),sb(pm  +  1,  k), . . . ,  sb((p  +  l)m  —  1,  k)]T . 
The  corresponding  measurement  equation  is  then 

u p(k)  =  s bp(k)  +  w p(k),  p  =  0, 1, . . . ,  P  -  1.  (3) 
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In  the  case  where  the  target  return  is  present  in  the  received  signal  the  state  space 
representation  is  defined  as: 


Sp(fe) 

=  4(k  -  1) 

(4) 

■£(*) 

=  sp(fc  -  1)  +  bp(fc) 

(5) 

b p{k) 

=  bp{k  -  1)  +  v2p(fc) 

(6) 

for  p  =  0, 1, . . . ,  P  —  1.  The  state  is  augmented  with  the  vector  s p(k)  =  [s*(pm,  k),st(pm  + 
1,  k), . . . ,  sl((p  4-  l)m  —  1,  k)]T  and  a  random  bias  bp(k)  which  is  used  to  account  for  the 
changes  corresponding  to  the  target  signal.  The  background  is  kept  constant  and  equal 
to  the  value  estimated  prior  to  target  detection  (see  Eq.  (3)).  For  the  above  augmented 
state  the  measurement  equation  is  defined  as: 

u p(k)  =  s bp{k)  +  s p(fc)  +  Wp(fc),  p  =  0, 1, . . . ,  P  -  1.  (7) 

We  note  that  vip(k)  in  (2)  and  V2P{k)  in  (6)  are  the  process  noise  vectors  and  wp(k)  in  (3) 
and  (7)  is  the  measurement  noise  vector,  for  p  =  0, 1, . . . ,  P  —  1.  They  are  all  m-dimensional 
vectors  whose  components  are  independent  identically  distributed  (iid)  Gaussian  random 
variables.  Also,  in  the  background-only  regime  the  target  signal  s*(n  =  0, . . . ,  N  —  1,  k)  is 
zero  so  that 

sj(/c)  =  [0]mxi  for  p  =  0,l,...,P-  1.  (8) 

Based  on  the  equations  (l)-(2)  and  (3)-(6),  a  set  of  P  Kalman  filters  is  used  to  com¬ 
pute  the  background  and  the  target  signal  estimates  from  the  noisy  observations.  In  the 
quiescent  phase  the  p-th  state  vector  xip  =  [sp]  is  recursively  estimated  as 

xi p(k)  =  Aixi p(k  -  1)  +  kip(A;)(up(fc)  -  HiAixip(A;  -  1))  (9) 

where  Ai  =  I  and  Hi  =  I,  and  I  is  an  identity  matrix  with  m  elements  on  the  main 
diagonal.  The  update  equations  in  this  case,  for  p  =  0, 1, . . . ,  P  —  1,  are 

Plp(fc|fc-1)  =  AiPip(fc  -  l)Af  4-  Qlp(fc)  (10) 

kip(fc)  =  Plp(fc|fc-l)Hf(HiPip(fc|A:-l)Hr  +  Rp)-1  (11) 

Plp(fc)  =  (I  —  klp(A:)H1)Pip(A:|A:  —  1)  (12) 

where  kip(A;)  is  a  Kalman  gain  vector,  Pip(k\k  -  1)  is  an  a  priory  error  covariance  ma¬ 

trix  with  the  dimensions  m  x  m,  Pip(&)  is  the  updated  error  covariance  matrix  with 
the  same  dimensions,  Qi p(k)  is  the  process  noise  covariance  matrix  and  Rp  =  er^I  is 
the  measurement  noise  covariance  matrix.  The  error  covariance  matrix  is  initialized  to 
Plp(0)  =  [0]mxm  and  the  process  noise  covariance  matrix  to  Qip(0)  =  (<r°lp)2I,  where 
(cr°lp)2  is  the  initial  estimate  of  the  variance  of  the  process  noise  sequence  in  (l). 

When  the  target  signal  is  present,  the  augmented  state  vector  is 

*2p  =  [  [sjf,  [s‘)T,  [bP]T]T  (13) 

and  its  updated  estimate  is 

X2p(fc)  =  A2x2p(fc  —  1)  +  k2p(fc)(up(fc)  -  H2A2x2p(fc  -  1))  (14) 
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where 

10  0 
A2=  Oil 
0  0  1 

Similarly  as  above,  gain  and  the  covariance  matrix  updates  for  p  =  0,1 ,P  —  1,  are 
obtained  as 


P2p(*|fc-1)  =  A2P2p(fc  -  1)A%  +  Q2p  (16) 

k :2p(k)  =  P2p(A;|A:-l)Hi,(H2P2p(fc|fc-l)Hi’  +  Rpr1  (17) 

P2p(fc)  =  (I  —  k2p(fc)H2)P2p(fc|fc  —  1).  (18) 

The  process  noise  covariance  matrix  Q2p  is  a  3 m  x  3m  matrix  defined  as 

0  0  0 

Q2p  =  0  0  0  .  (19) 

[  0  0  "Ip1  . 


The  augmented  state  is  initialized  at  the  spatial  position  k  —  ko  by  setting  the  bias  to 
zero  and  the  error  covariance  matrix  P2p(/co)  to 

‘  PiP(fco  -  1)  0  O' 

P2p(fco)=  0  0  0.  (20) 

0  0  0 

3  Target  Detection 

While  operating  in  the  quiescent  state  the  algorithm  uses  a  x2  test  based  on  the 
measurement  prediction  error  or  innovation  up(k ),  to  detect  the  position  of  a  possible 
target.  In  particular,  it  is  expected  that  the  presence  of  the  target  return  manifests  itself 
as  a  “large”  innovation.  Under  the  hypothesis  Ho  it  is  assumed  that  the  GPR  trace 
u(n  =  0, . . . ,  N  -  1,  k)  contains  only  the  ground  return  sb(n  =  0, . . . ,  N  -  1,  k)  plus  noise. 
The  alternative  hypothesis  H\  is  that  the  target  signal  st(n  —  0, . . .  ,N  —  1,  k)  is  also 
present.  For  each  of  the  P  trace  sections  of  length  m,  p  —  0, 1, . . . ,  P—  1,  the  measurement 
prediction  error  vp(k)  is  computed  as 

vp(k)  =  u p{k)  -  up{k\k  -  1)  (21) 

where 

up{k\k  -1)  =  HiAiXi p(k  -  1)  (22) 

is  the  measurement  prediction.  The  measurement  prediction  covariance  matrix  is  updated 
as 

S p{k)  =  HjPi p(k\k  -  l)Hf  +  Rp  (23) 

and  used  to  compute  the  normalized  innovation  squared  (NIS) 

ep(k)  =  vp{k)TSp{k)-1up{k).  (24) 
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Under  the  hypothesis  H0  the  variables  ep(k),  p  =  0, 1, . . . ,  P  —  1,  are  x2  distributed  with 
m  degrees  of  freedom.  Each  of  the  P  trace  sections  is,  then,  tested  separately  using  the 
X 2  test  with  the  significance  level  a,  so  that 

Pr(ep(k)  >  x2m,a)  =  a-  (25) 

In  the  case  when  Hq  is  rejected  for  at  least  Kq  of  the  total  of  P  sections,  the  hypothesis  Hq 
for  the  whole  trace  is  rejected  and  the  alternative  hypothesis  Hi  is  accepted.  To  make  the 
final  decision  that  the  target  is  present,  and  to  change  to  the  augmented  state  model,  the 
hypothesis  Ho  needs  to  be  rejected  for  at  least  K\  consecutively  measured  GPR  traces. 
After  the  target  detection  is  declared  at  the  spatial  position  k,  the  spatial  position  at 
which  the  augmented  state  is  initialized,  ko,  and  which  is  assumed  to  correspond  to  the 
target  onset,  is  determined  as  ko  =  k  —  K\  —  KT.  Here,  Kq ,  K\  and  KT  are  suitably  chosen 
constants.  If  KT  =  0,  ko  equals  to  the  spatial  position  where  detection  first  occurred.  If 
we  set  Kt  >  0  we  assume  that  the  target  started  before  it  was  first  detected. 

It  is  desirable  that  the  augmented  state  model  stays  active  until  the  end  of  the  section 
of  GPR  traces  that  contains  the  target  return  is  detected,  when  the  algorithm  should 
switch  back  to  the  quiescent  state  model.  However,  the  above  target  detection  technique 
determines  the  beginning  of  the  target  section  in  the  data  and  a  similar  approach  cannot 
be  devised  to  define  the  end  of  this  segment  as  easily.  One  solution  is  to  run  the  proposed 
algorithm  over  the  target  in  two  opposite  directions  in  order  to  define  the  position  and  the 
width  of  the  target  segment.  For  simplicity,  in  this  report  we  apply  the  algorithm  in  one 
direction  only  and  assume  that  the  width  of  the  target  segment  is  constant  and  known  in 
advance.  This  assumption  is  reasonable  when  the  size  of  the  target  and  the  velocity  of  the 
antenna  are  known. 

Though  here  we  used  NIS  for  the  detection  statistic,  we  note  that  the  detection  statistic 
can  also  be  computed  as  a  moving  sum  of  NIS  over  the  sliding  window  of  length  s, 

ep(fc)  =  (26) 

j=k—s+ 1 

The  variable  ep(k)  is  x2  distributed  with  sm  degrees  of  freedom.  The  detection  statistic 
defined  as  in  (26)  is  suitable  for  noisier  backgrounds. 


4  Background  Adaptation 

The  GPR  background  signal  is  measured  as  the  average  over  a  volume  of  the  ground 
and  its  statistics  can  be  assumed  to  have  slow  spatial  variation.  To  enable  our  algorithm 
to  adapt  to  such  variations  we  use  an  approach  motivated  by  the  continuous  noise  level 
(CNL)  adjustment  technique  presented  in  Bar-Shalom  and  Li  [1].  In  particular,  the  process 
noise  covariance  matrix,  Qip(k),  that  characterizes  the  background  is  finely  tuned  at  each 
spatial  position  k  by  multiplying  it  by  a  scaling  factor  (f>p{k), 

Qip(k)  =  <t>p{k)Qip{k  -  1).  (27) 

This  is  done  for  each  p  =  0, 1, . . . ,  P  —  1  separately.  At  each  spatial  position  k  the  value 
of  the  scaling  factor  cf)p(k)  is  determined  by  comparing  NIS  ep{k)  to  a  set  of  thresholds. 
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These  thresholds,  denoted  as  ei,  e2 L  and  e2w ,  are  defined  so  as  to  satisfy  the  following 
equations 


Pr(ep(k)  <  ei) 

=  i-m 

(28) 

Pr{ep(k)  <  e2J 

=  1  -  V2L 

(29) 

Pr(ep(k)  <  e2ll) 

=  i-V2„ 

(30) 

where  771  >  r}2L  >  r)2H  >  a  and  Pr(-)  is  x2  distribution  with  m  degrees  of  freedom.  Ideally, 
under  the  background-only  assumption,  we  wish  to  define  the  process  noise  covariance 
matrix  Q\p{k)  so  that  <  ep(k)  <  e2l,  in  which  case  the  scaling  factor  is  set  to  <f>p(k)  =  1. 
If  ep(k)  <  €1,  the  process  noise  level  is  too  high  and  needs  to  be  lowered.  We  then  set 
(f>p(k)  =  <pu  where  <f>  1  <  1.  Otherwise,  if  e2i  <  ep(k)  <  e2/f,  the  process  noise  covariance 
matrix  is  increased  by  setting  the  scaling  factor  to  4>p(k )  =  4> 2,  where  </> 2  >  1-  In  the  case 
when  €2h  <  Cp(fc)  <  Xm,a  we  are  not  sure  weather  our  background-only  assumption  is 
still  valid  so  the  process  noise  level  is  left  unaltered,  and  4>p(k)  =  1.  The  same  is  done  for 
tp{k)  >  Xm,a •  The  constants  (f>\  and  </>2  are  chosen  to  be  close  to  1  and,  i.e.,  <f>\  =  0.98 
and  (f>2  =  1-02.  The  decision  procedure  for  setting  the  value  of  the  scaling  factors  4>p(k)  at 
each  spatial  position  k  for  which  the  algorithm  operates  using  the  quiescent  state  model 
and  for  p  —  0, 1, . . . ,  P  —  1,  is  summarized  as  follows  (see  also  Figure  1): 

step  1:  set  p  —  0 

step  2:  if  ep(k)  <  ei  set  <j>p{k)  =  cj)\  goto  step  3 

if  C2L  <  eP(k)  <  e2li  set  <pp(k)  =  fc  goto  step  3 
if  d  <  tp(fc)  <  e2i  or  ep(k)  >  xh,a  set  <f>p{k)  =  1 
step  3:  if  p  <  P  set  p  =  p  +  1  goto  step  2 

5  Simulation  Results  and  Discussion 

The  Kalman  filter-based  technique  described  in  Sections  2-4  has  been  tested  using 
several  minelike  targets  buried  in  different  soils.  The  targets  were  surrogate  anti-personnel 
landmines,  modelled  after  a  number  of  anti-personnel  blast  mines  with  non-metallic  casings 
[7],  and  poly-vinyl  chloride  (PVC)  and  stainless  steel  cylinders  of  different  sizes.  The 
surrogate  landmines  are  made  of  a  plastic  pipe  filled  with  paraffin  wax  and  small  metal 
parts.  The  data  was  collected  by  means  of  an  FR-127-MCSB  impulse  GPR  developed  by 
CSIRO  [6].  The  system  used  a  bistatic  bow-tie  antenna  to  transmit  pulses  with  1.4  GHz 
center  frequency  and  the  bandwidth  of  1.0  GHz,  and  collected  127  soundings  composed 
of  512  samples  of  12  bit  accuracy  per  second.  The  length  of  each  trace  was  12  ns.  To 
compensate  for  the  ground  attenuation  and  signal  spreading  the  data  was  taken  using  the 
dynamic  gain  slope  of  20  dB.  In  the  process  of  measurement  the  antenna  was  suspended 
from  a  track  along  which  it  was  driven  by  a  stepper  motor  at  constant  velocity.  More 
information  about  the  measurement  process  and  conditions  can  be  found  in  [3]. 

The  proposed  technique  worked  well  and  was  found  promising  for  detecting  weak  tar¬ 
gets.  One  such  example  is  shown  in  Figure  2.  Figure  2  (a)  shows  the  original  GPR 
plot  that  contains  two  targets,  the  surrogate  landmine  ST-AP(l)  [7]  and  a  PVC  cylinder, 
buried  in  dry  sand  at  the  depth  of  approximately  1  cm.  The  dimensions  of  the  targets 


(31) 


(32) 
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axe:  ST-AP(l):  diameter  5.2  cm  and  height  4.2  cm;  PVC  cylinder:  diameter  of  10  cm  and 
height  5  cm.  The  target-return  component  of  the  GPR  signal  estimated  by  the  proposed 
algorithm  is  shown  in  Figure  2  (b). 

Since  the  individual  GPR  traces  were  oversampled  in  the  process  of  measurement, 
each  of  them  was  decimated  by  factor  2  prior  to  the  processing.  The  parameters  of  the 
algorithm  used  to  obtain  the  results  in  Figure  2  were:  the  width  of  the  horizontal  strips 
of  GPR  data  processed  by  one  Kalman  filter  m  =  32  and  the  number  of  such  horizontal 
strips  P  =  8.  The  significance  level  for  the  x2  test  was  a  =  10-°.  Other  parameters  used 
in  the  detection  procedure  were  Kq  =  1,  K\  =  5  and  Kr  =  5,  and  for  the  background 
adaptation  r)\  =  0.99,  r]2L  =  0.4,  rj2H  =  5a.  For  all  P  =  8  trace  sections  the  variances  of  the 
components  of  the  measurement  noise  vectors  wp  were  identical  and  equal  to  aw  which  was 
estimated  off  line,  similarly  as  the  variances  of  the  components  of  the  process  noise  vectors 
in  (6).  Since  the  intention  was  to  detect  shallowly  buried  targets,  the  target  detection 
using  x2  hypothesis  testing  was  performed  only  for  the  top  6  horizontal  data  strips  where 
the  specular  return  from  the  targets  was  expected  to  be  maximum.  The  more  accurate 
detection  results  were  obtained  for  the  larger  values  of  the  parameter  m.  The  drawback 
was  that  the  computational  complexity  increased  with  the  size  of  the  matrices  used  in 
the  algorithm.  In  our  experiments  m  =  32  has  been  found  to  enable  very  accurate  target 
detection  with  the  acceptable  computational  complexity.  The  program  is  implemented 
using  Matlab  Numeric  Computation  and  Visualisation  Software  and  is  listed  in  Appendix 
A.  The  flow  chart  of  the  algorithm  is  shown  in  Figure  3. 

Figure  4  (a)  shows  a  row  from  the  data  in  Figure  2  (a),  which  represents  the  signal 
recorded  at  one  particular  depth  and  at  consecutive  spatial  positions  (solid  line).  The 
corresponding  background  signal  estimated  using  the  Kalman  filter-based  approach  is  also 
shown  (dashed  line).  Figure  4  (b)  shows  the  estimated  target  signal.  The  augmented  state 
model  was  activated  over  the  sections  of  the  signal  between  the  vertical  lines,  as  indicated 
in  the  figure.  It  can  be  seen  that  the  background  follows  the  received  signal  at  the  spatial 
positions  where  the  target  signal  is  zero,  that  is,  when  the  quiescent  state  model  is  used. 
At  the  spatial  positions  where  the  target  signal  is  present,  i.e.,  when  the  augmented  state 
model  is  active,  the  background  signal  stays  at  the  level  estimated  before  the  detection. 

For  comparison,  Figure  5  shows  the  data  in  Figure  2  (a)  processed  by  background 
subtraction,  where  the  background  signal  b(n)  is  estimated  as  the  mean  of  the  entire  data 
ensemble  in  Figure  2  (a)  (i.e.,  b(n)  =  IC^io1  ui fc(n),  n  =  0, 1, . . . ,  N,  and  M  is  the  overall 
number  of  traces  in  the  ensemble)  and  subtracted  from  each  trace  in  the  ensemble.  In 
Figure  5  it  can  be  seen  that,  using  this  approach,  not  all  clutter  has  been  removed  from 
the  data. 


6  Conclusions  and  Future  Work 


The  signatures  of  shallowly  buried  targets  measured  using  Ground  Penetrating  Radar 
are  usually  obscured  by  the  return  from  the  air-soil  interface.  In  this  report  a  Kalman 
filter-based  approach  has  been  used  to  obtain  the  estimate  of  this  background  signal  and 
to  separate  it  from  the  target  return.  In  the  absence  of  the  target  the  filter  operates 
using  a  quiescent  state  model  in  which  it  computes  the  background  estimate.  When  the 
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target  is  detected,  new  state  components  are  added  to  allow  for  the  change  of  the  signal 
corresponding  to  the  target  return.  The  augmented  state  model  is  used  until  it  reverts  to 
the  quiescent  model  by  another  decision. 

The  magnitude  of  the  background  signal  at  one  particular  depth,  and  for  incrementally 
updated  spatial  positions,  was  modelled  as  a  random  walk.  The  noisy  measurements  of 
the  background  signal  were  represented  using  a  state  space  form  and  a  set  of  Kalman 
filters  was  used  to  obtain  the  recursive  state  estimation.  At  the  spatial  positions  where  a 
target  was  detected  the  new  state  component  modelled  as  the  stochastic  bias  was  added. 
This  component  was  used  to  compute  the  estimate  of  the  unknown  target  signal. 

To  detect  the  changes  in  the  signal  that  are  related  to  the  presence  of  the  target  a 
statistic  based  on  the  measurement  innovation  is  used  in  conjunction  with  an  appropriate 
statistical  hypothesis  testing  procedure. 

The  proposed  algorithm  was  tested  using  the  data  containing  several  minelike  targets 
and  have  shown  promise  for  processing  GPR  data  with  weak  shallow  target  responses. 

The  future  work  will  involve  using  a  more  elaborate  Multiple  Model  (MM)  Kalman 
filter  approach  [1].  The  assumption  of  the  MM  method  is  that  the  signal  obeys  either  of  the 
two  possible  models  (i.e.,  either  the  target  return  is  present  in  the  signal  or  not,  similarly 
as  described  in  this  report).  It  computes  the  probability  that  each  of  the  models  is  active  at 
one  particular  spatial  position  using  a  Bayesian  approach  and  can  incorporate  Markovian 
assumption.  It  is  expected  that  more  accurate  estimation  of  both  the  background  and 
the  target  signal  can  be  obtained  using  this  method.  It  will  also  include  extending  the 
approach  presented  in  [4]. 
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Figure  1:  Thresholds  e\,  eiL,  £2H  and  Xm,a  and  the  corresponding  settings  of  the  scaling 
factors  (f>p(k).  Function  fx^  is  a  x2  probability  density  function. 
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Figure  2:  (a)  The  original  GPR  data  with  the  surrogate  landmine  ST-AP(l)  (left)  and 
the  PVC  cylinder  (right)  buried  in  dry  sand  at  the  depth  of  1  cm;  (b)  the  target-return 
component  of  the  GPR  signal  estimated  by  the  proposed  algorithm.  The  positions  of  the 
targets  are  indicated. 
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Figure  3:  Flow  chart  of  the  algorithm. 
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(a) 


(b) 

Figure  4'  (a)  A  row  from  the  data  plot  in  Figure  2  (a)  which  represents  the  signal  recorded 
at  one  depth  and  at  consecutive  spatial  positions  (solid  line)  and  the  corresponding  back¬ 
ground  signal  estimated  using  the  proposed  algorithm  (dashed  line);  (b)  the  estimated  target 
signal.  The  augmented  state  model  was  active  at  spatial  positions  as  indicated. 
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Figure  5:  Background-subtracted,  data  in  Figure  2  (a),  where  the  background  signal  is 
estimated  as  the  mean  of  the  entire  ensemble  of  GPR  signals  and  subtracted  from  each 
signal  in  the  ensemble. 
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i***********^***** ************************************************ ******** 

*/. 

X  Function:  [M,Ml,Bl]=smooth_rw8(D,dim)  ; 

*/, 

*/,  This  function  applies  Kalman  filtering  based  on  the  variable 
*/,  state  dimension  method  as  expalined  in  thie  report. 

7. 

D  two-dimensional  plot  of  GPR  data 

dim  -  width  of  data  strips  processed  by  one  Kaman  filter 


'/,  Input : 

% 

y. 

*/,  Output: 

y. 
y. 
y. 

'/,  Author: 

y. 
x 

X  Date: 

X 

**])c:|c:|c******************************************************************** 


M  -  estimated  background  signal 
Ml  -  estimated  target  signal 
B1  -  estimated  bias 

Dragana  Carevic 

Surveillance  Systems  Division,  DSTO 
18/11/1998 


function  [M,M1 ,Bl]=smooth_rw8(D,dim) ; 


X  Filter  out  the  low-pass  component  of  the  signal 
qmf=MakeONFilter(’ Coif let’ ,1) ; 

[X,Y]=size(D) ;N=4;k=l;  step=4; 

Xdisp([’  Step  =  ’ ,num2str(step)3 ) ; 

for  i  =  1 :step:Y-step+l 
d=D(: ,i) ; 

d=decimate(d,2, ’fir') ; 

*/,d  =  DownDyadLo(d’ ,qmf )  ’ ; 
wave_coef  =  FWT_PO(d,N,qmf ) ; 

wave_coef  (1:(2"N))  =  zeros (size (wave_coef  (1: (2~N)))) ; 
dfilt( : ,k)=  IWT_PO(wave_coef ,N,qmf ) ; 
k=k+l ; 
end; 
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numChsqAverage  =  1; 

chisqPrev  =  zeros (1 .numChsqAverage) ; 

anythingDetected  =  0;  numLess  =  0; 

Detection  =  0;  Thres  =  -5; 
detected=0;  detectedRequired=5 ; 

LookAhead  =  60;  7,  distance  to  be  looked  at  in  detection 
Num  =  ceil (LookAhead/ step)  -  1; 

NumDetectedTraces  =  0;  NumBkgTraces  =  0; 

PrevDetection  =  1; 

v=40*step;  7.  std  of  the  measurement  noise  (cov  matrix  R) 

w0=20*step;  7.  std  of  the  process  noise  (for  the  background)  (cov  matrix  Q) 

[XF,YF]=size(dfilt) ; 

m  =  f  ix(XF/dim) ;  7.  number  of  strips  to  be  processed 
Z  =  zeros  (dim,  dim) ;  7.  zero  matrix  of  size  dim  x  dim 
Xi  =  eye  (dim);  7.  identity  matrix  of  size  dim  x  dim 
H  =  eye (dim) ; 

A  =  eye (dim) ; 

R  =  v.  ~2*eye(dim)  ;  7,  measurement  covariance  matrix  dim*dim 


for  i  =  l:m 
Q(i,:,0  = 

P(i,:,0  = 

K(i, : , :)  = 
end; 


w0.~2*eye(dim) ;  7.  process  covariance  matrix  dim*dim 

zeros  (dim,  dim)  ;7.eye  (dim)*5000“2 ; 

7.  initialize  error  covariance  matrices  for  m  strips 
7.  of  dimensions  m  x  dim  x  dim 
zeros (dim, dim) ; 

7,  initialize  gain  matrix  for  m  strips  m  x  dim  x  dim 


M( : ,1)  =  dfilt( : ,1) ; 

7,  initialize  vectors  M(:,i)  to  values  defined  in  dfilt(:,l) 
7,  M  corresponds  to  the  background 
M1(1:XF,1)  =  zeros(XF.l); 

7,  initialize  vectors  Ml(:,i)  to  zero 
l  Ml  corresponds  to  the  signal 
B1(1:XF,1)  =  zeros(XF.l); 

7,  initialize  vector  Bl(:,l)  to  zero. 

7.  B1  corresponds  to  the  signal 

j  =  2; 

while  j  <  YF-1 

PrevDetection!nThisTrace=0; 
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disp([’  Trace  =  ’  ,num2str(j)]  )  ; 

disp([’  NumDetectedTraces  =  ’ ,num2str(NumDetectedTraces)] ) ; 
disp([’  NumBkgTraces  =  ’ ,num2str(NumBkgTraces)] ) ; 
disp([’  PrevDetection  =  ’ ,num2str(PrevDetection)] ) ; 

disp([’  PrevDetectionlnThisTrace  =  ’ ,num2str(PrevDetectionInThisTrace)] ) ; 
if  Detection  ==  0 

NumBkgTraces  =  NumBkgTraces  +1; 
if  NumBkgTraces  >  200/step 
PrevDetection=0 ; 
end; 


for  i  =  l:m  */,  for  the  m-th  strip  at  column  j 

s=[M(dim*(i-l)+l :dim*i, j-1) *] ’ ; 

'/,  this  is  previous  state  vector 
7,  based  on  which  we  are  making 
*/,  current  state  prediction 


7,  PREDICTION 

sp  =  A*s;  7.  this  is  current  state  prediction 
PP(: , :)-P(i, : , ;  QQ  =  Q(i,:,:>; 

PP=A*PP*A’+  QQ;  7.  prediction  error  covariance  matrix 
mp  =  H*sp;  7.  measurement  prediction 

S  =  H*PP*H’+R;  7.  measurement  prediction  covariance  matrix 

7.  CHI  SQUARED  TEST  FOR  TARGET  DETECTION 
for  1  =  1:1 

error  =  dfilt(dim*(i-l)+l:dim*i,j+l)  -  mp; 

7.  measurement  prediction  error 


chisq  =  error’ *inv(S)*error; 
for  wtq=l :numChsqAverage-l 
chisqPrev(wtq)=chisqPrev(wtq+l) ; 
end; 

chisqPrev(numChsqAverage)=chisq; 
chisq=sum(chisqPrev) ; 

7.  chi-square  variable  of  dim  degrees  of  freedom 
7.  mean  previously  computed  numChsqAverage  values 

Prob(i, j)=l-chi2cdf (chisq, dim*numChsqAver age)  ; 

7.disp( [’Meras .  innovation  :  Prob  =  *  ,num2str((Prob(i,  j)))] ) ; 

if  j>f ix(200/step)  ft  (Prob(i,j))  <=  10~Thres  ft  i<  m/2  ft  ... 
PrevDetection  ==  0  ft  PrevDetectionlnThisTrace  ==  0 
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disp(C’ DETECTION  at  strip  =  ’ ,num2str(i)] ) ; 
disp([’Prob  =  ’ ,num2str((Prob(i, j)))] ) ; 

if  detected  "=  0  &  detectedTrace (detected)  ==  j 
break; 
end; 

anythingDetected  =  1 ; 

PrevDetectionlnThisTrace  =  1; 
if  detected==0 

%  if  this  is  first  detected  trace 
%  initialize  new  detection 

detectedTrace (l)=j ; 
detected  =  1; 

else  '/,  check  if  there  are  previous  consecutive  traces 

if  detectedTrace (detected) +1  ~=  j 

%  if  not,  reset  variable  detected  to  zero 

detected  =  1; 

detectedTrace (l)=j ;  %  initialize  new  detection 

else  7,  if  the  traces  are  consequtive  increment 
7,  variable  detected  and  check  if  there  is 
7.  detectedRequired  number  of  them 

detected  =  detected  +i; 
detectedTrace (detected) =j ; 

end; 

end; 

if  detected  ==  detectedRequired 
Detection  =1; 

DetectionTrace=j-l ; 

j  = j - (detectedRequired+9+2*numChsqAverage) ; 

end; 

disp( [’detected  =  ’ ,num2str (detected)] ) ; 

7,  If  detection  occur  switch  to  other  (compund)  model 
7,  In  particular  break  out  of  this  loop  and  use  j-l-th 
7,  estimate  as  the  initial  value  for  the  new  model 

break; 
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end; 

end; 

if  Detection  ==  1  break;  end; 

7.  MEASUREMENT  UPDATE 
K  =  PP*H’/ (H*PP*H’+R) ; 

s  =  sp  +  K*(dfilt(dim*(i-l)+l:dim*i, j)  -  H*sp) ; 

7,  update  state  vector 
P(i, : , : )  =  PP  -  K*H*PP; 

7*  update  error  covariance  matrix  for  i-th  strip 
M(dim*(i-1)+1 :dim*i, j)=s(l:dim) ; 

7.  move  current  m  vector  to  matrix  M(:,i) 
Ml(dim*(i-l)+l:dim*i, j)=zeros(dim, 1) ; 

Bl(dim*(i-1)+1 :dim*i, j)=zeros(dim,l) ; 

7.  since  no  detection  occured,  B1  and  Ml  are  zero 

7.  BACKGROUND  ADAPTATION 

if  min(Prob(i, j))  >  0.99  ft  anythingDetected  ==  0 
QQ  =  0.98  *QQ; 
end; 

if  min(Prob(i,j))  <  0.4  ft  anythingDetected  ==  0  ft  ... 
min(Prob(i , j))  >  5*10~Thres 

numLess  =  numLess  +1; 
if  numLess  >=  1 

QQ  =  1.02  *  QQ; 

end; 

else  numLess  =0; 

Q(i> : > : )=QQ ; 
end; 

anythingDetected  =  0; 
if  Detection  ==  0 

j=j+l  ; 

end; 

end  '/,  if  Detection  ==  0 
if  Detection  ==  1 

if  PrevDetection  ==  0 

7.  Initialisation  of  parameters  for  compond  model 
wl  =  (70/5)*step; 
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III  =  eye(3*dim); 

I  =  eye (dim) ; 

HI  =  [I  I  Z]  ; 

Al  =  [I  Z  Z;  Z  I  I;  Z  Z  I]; 

Q1  =  [Z  Z  Z;Z  Z  Z;  Z  Z  vl~2*I]  ; 

'/,  Initialise  new  values  for  PI  and  K1 
7.  Use  old  values  P  and  K  for  the  background 
for  i  =  1  :m 
PP(:,:)»P(i,:» :); 

Pl(i, : , =  [PP  Z  Z;Z  Z  Z;Z  Z  Z] ; 
end; 

PrevDetection  =  1 ; 
end;  7.  if  PrevDetection  ==  0 
NumDetectedTraces  =  NumDetectedTraces  +1; 


for  i  =  l:m 


si  =  [M(dim*(i-1)+1 :dim*i, j-1) ’  Ml(dim*(i-l)+l:dim*i, j-i) ’  ... 

Bl(dim*(i-1)+1 :dim*i, j-1) ’] ’ ; 


7.  PREDICTION  STEP 
spl  =  Al*sl; 

PPl(:,:)-Pl(i,:.:>; 

PP1  =  A1*PP1*A1’+Q1; 

7.  MEASUREMENT  UPDATE 

K1  =  PP1*H1 ’ / (H1*PP1*H1 ’+R) ; 

si  =  spl  +  Kl*(dfilt(dim*(i-1)+1 :dim*i, j)  -  Hl*spl) ; 

7.  update  state  vector 
Pl(i, : , =  PP1  -  K1*H1*PP1; 

7,  update  error  covariance  matrix  for  the  i-th  strip 
M(dim*(i-l)+l:dim*i, j)  =  sl(l:dim); 

Ml(dim*(i-l)+l:dim*i, j)  =  sl(dim+l:2*dim) ; 

Bl(dim*(i-1)+1 :dim*i, j)  =  si (2*dim+l :3*dim) ; 

end;  7.  for  i=l:m 

Prob(: ,j)=Probl(: ,j) ; 

j=j+i; 


if  NumDetectedTraces  >  300/step 
detected  =  0; 

Detection  =  0; 
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NumBkgTraces  =  0; 
NumDetectedTraces  =0; 

end; 

end;  */,  if  Detection  ==  1 
end; 
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